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ABSTRACT 

Context. We study turbulent convection during the core helium flash close to its peak by comparing the results of two and three- 
dimensional hydrodynamic simulations. 

Aims. In a previous study we found that the temporal evolution and the properties of the convection inferred from two-dimensional 
hydrodynamic studies are similar to those predicted by quasi-hydrostatic stellar evolutionary calculations. However, as vorticity is con- 
served in axisymmetric flows, two-dimensional simulations of convection are characterized by incorrect dominant spatial scales and 
exaggerated velocities. Here, we present three-dimensional simulations that eliminate the restrictions and flaws of two-dimensional 
models, and that provide a geometrically unbiased insight into the hydrodynamics of the core helium flash. In particular, we study 
whether the assumptions and predictions of stellar evolutionary calculations based on the mixing-length theory can be confirmed by 
hydrodynamic simulations. 

Methods. We use a multidimensional Eulerian hydrodynamics code based on state-of-the-art numerical techniques to simulate the 
evolution of the helium core of a 1.25M Q Pop I star. 

Results. Our three-dimensional hydrodynamic simulations of the evolution of a star during the peak of the core helium flash do 
not show any explosive behavior. The convective flow patterns developing in the three-dimensional models are structurally different 
from those of the corresponding two-dimensional models, and the typical convective velocities are smaller than those found in their 
two-dimensional counterparts. Three-dimensional models also tend to agree better with the predictions of mixing length theory. Our 
hydrodynamic simulations show the presence of turbulent entrainment that results in a growth of the convection zone on a dynamic 
time scale. Contrary to mixing length theory, the outer part of the convection zone is characterized by a sub-adiabatic temperature 
gradient. 
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1. Introduction 

The core helium flash is the most violent event in the life of 
a star with an initial mass between approximately 0.7M Q and 
2.2M Q (Sweig art & Gross|1978| l. Electron-degeneracy in the he- 
lium core at the time of the flash leads to a thermonuclear run- 
away producing, at its peak, large amounts of energy (~ 10 10 L G ) 
within the stellar core over a very short period of time (~ days). 
The temperature rises up to several 10 8 Kelvins until the degen- 
eracy of the electron gas is eventually lifted. Nevertheless, the 
event seems not to be catastrophic for the star. It results only in a 
slow expansion of the helium core (typically with a few m s _1 ), 
since energy transport by turbulent convection, heat conduction, 
and radiation seems to be able to deliver most of the flash energy 
quiescently from the stellar interior to the outer stellar layers. 
However, computations of the flash have a confusing history pre- 
dicting either severe explosions (Edwards 1969; Zimmermann 
l^[Wiere|1976||Wickett|1977ljCole & Deupree|1980)M981 



Deupree & Cole||1983| |Deupree||1984a|b| |Cole et aL|[1985 



or a quiescent behavior ( |Deupree||1996| |Dearborn et al.| 2006 
|Lattanzio et al.|2006) . 

Previous two-dimensional hydrodynamic simulations of tur- 
bulent convection within the helium core during the peak of 
the core helium flash (Mocak et al. 2008) support a quies- 
cent scenario in agreement with stellar evolutionary calculations. 
However, they also showed that the convection zone, powered by 



helium burning, is characterized by convective velocities that are 
roughly four times larger than those predicted by mixing-length 
theory, and that the width of the convection zone grows on a dy- 
namical timescale. This may lead to mixing of hydrogen into the 
helium core as it is known from one-dimensional simulations 



of extremely metal-poor stars (Fujimoto et al. 1990 Schlattl 



et al.|2001[|Cassisi et al.|2003||Weiss et al.||2004[ |Campbell & 



Lattanzio[2 008) 



It is well known that two-dimensional hydrodynamic sim- 
ulations of turbulence are seriously biased due to the imposed 
symmetry restrictions. Opposite to 3D flows, the turbulent ki- 
netic energy increases from small to large scales in 2D simula- 
tions, i.e., the energy cascade to smaller length scales character- 
istic of turbulent flows is not reproduced (Canuto 2000 Bazan & 
Arnett|1998) . Hence, the mean convective velocities, the amount 
of overshooting, and the size of turbulent structures are too large. 
Thus, as pointed out already by e.g., |Muthsam et aL ( 1995 1 and 
Baza n & Arnett| ([1998 ), three-dimensional simulations are re- 
quired to validate the predictions of two-dimensional simula- 
tions. 

With increasing computational capabilities, the importance 
of multidimensional hydrodynamic simulations in stellar evo- 
lution studies grows rapidly because they are essentially "pa- 
rameter free". In contrast, canonical one-dimensional stellar 
evolutionary calculations involve free parameters like the mix- 
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Table 1. Some properties of the initial model: total mass M, stellar population, metal content Z, mass Mhc and radius Rn e of the 
helium core (X( 4 He) > 0.98), nuclear energy production in the helium core Lfj e , temperature maximum T max , radius r max , and 
density p max at the temperature maximum. 
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Fig. 1. Left panel: Temperature (in units of 10 7 K) distribution of the mapped (dashed) and stabilized (solid) initial model displayed 
together with the density (in units of 10 5 g cm* 3 ) and entropy (in k^ per baryon) stratification of the stellar evolution model (dotted 
and dash-dotted), respectively. Right panel: Chemical composition of the initial model. The dotted vertical lines mark the edges of 
the convection zone. 



ing length or the overshooting distance (Bohm-Vitense 1958; 



Canuto & Mazzitelli 1991 1, which are chosen in an appropri 
ate manner to fit observational data (Montalba n et al.| |2004). 
Comparison of the results obtained with both approaches is cru- 
cial, because it allows one to constrain, validate or disprove the 
free parameters used in stellar evolutionary calculations, and 
because it can provide a clue to the applicability limits of the 
canonical (ID) treatment (IBazan & Arnett 1998] IKercek et al. 



T998l[T999||Asida & Arnett|2000| |Herwig et al.|2006UMeakin 
& Arnett|20061|Dearborn et al.|2006[|Arnett et al.|2007t|Meakin| 
& Arnett|2007b| l. 

In the following we present an investigation of the core he- 
lium flash by means of two and three-dimensional hydrody- 
namic simulation using state-of-the-art numerical techniques, a 
detailed equation of state, and a time-dependent gravitational po- 
tential. Our hydrodynamic models are characterized by a con- 
vectively unstable layer (the convection zone) embedded in be- 
tween two stable layers composed of several chemical nuclear 
species and of a partially degenerate electron gas. Similar sys- 
tems were studied in the past by many authors e.g.,|Hurlburt et al.| 
T986) [1994]>; r 



and Brummel 



Muthsam et al. ( 1995 >; Singh et al. 



( |1995||1998| l 



scribe our hydrodynamics code in Sect. 3, and discuss and com- 
pare the results of the two and three-dimensional simulations in 
Sect. 4. In particular, we compare the temporal evolution, the ki- 
netic energy density, the power spectra, the thermodynamic fluc- 
tuation amplitudes, the stability of the flow structures, and the 
turbulent entrainment at the convective boundaries of our mod- 
els. We discuss our results in the context of the predictions of 
mixing-length theory in Sect. 5. Finally, we present a high res- 
olution 2D simulation covering almost 1 .5 days of core helium 
flash evolution in Sect. 6. We end with a summary of our findings 
and some conclusions in Sect. 7. 



2. Initial model 

The initial model was obtained with the stellar evolution code 
Garstec ( [Weiss & Schlattl||2000"l |2007] >. Some of its properties 
are listed in Table[T] and the distributions of temperature, density 
and composition are displayed in Fig.[T] The model possesses a 
white dwarf like degenerate structure with an off-center temper- 
ature maximum resulting from plasma and photo-neutrino cool- 



et al.| ( |2002| ) assuming, however, that the stellar ing. Its central density is 7 x 10 g cm . At the outer edge of the 

isothermal region in the center of the helium core the tempera- 
ture jumps up to T max ~ 1.7 x 10 8 K, and the adjacent convec- 
tion zone has a super-adiabatic temperature gradient. The core is 
predominantly composed of 4 He (mass fraction X( 4 He)> 0.98), 
and some small amounts of 'H, 3 He, 12 C, 13 C, 14 N, 15 N and 16 0, 
respectively . For our hydrodynamic simulations we adopt the 
abundances of 4 He, 12 C and ls O as the triple-a reaction dom- 
inates the nuclear energy production during the flash. The re- 



matter is composed of a single ideal Boltzmann gas. This gives 
extra relevance to our simulations because they allow us to study, 
e.g., the dependence of turbulent entrainment and the structure 
of convective boundary layers on the composition of the stel- 
lar gas, and on the composition gradients present in the stellar 
model. 

We introduce the stellar model used as input for our multi- 
dimensional simulations in the next section. We then briefly de- 
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maining composition is assumed to be adequately represented 
by a gas with a mean molecular weight equal to that of 20 Ne. 

As our multidimensional hydrodynamics code Herakles (see 
next section) utilizes an Eulerian computational grid, our initial 
data for the hydrodynamic simulations are obtained by polyno- 
mial interpolation from the stellar model which was computed 
on a Lagrangian grid. The initial hydrodynamics model obtained 
in this way is not in hydrostatic equilibrium, because the equa- 
tion of state of Timmes & Swesty ( 2000| l implemented in our 
hydrodynamics code differs from that of Rogers et al. (1996) 
implemented in the Garstec code (for a given density and tem- 
perature in the initial stellar model, the pressure differs by about 
1%). Stabilization of this initial hydrodynamics model resulted 
in a small decrease of the temperature compared to the stellar 
model (PigJTJ. 

3. Hydrodynamic code 

The hydrodynamic simulations were performed with an en- 
hanced version of the grid-based code Herakles (Mocak et al. 
2008). The adopted mathematical model implemented in the 
code consists of the Euler equations in spherical coordinates 
(r, 6, and <f>; see Appendix [A]) coupled to the source terms de- 
scribing thermal transport, self-gravity, and nuclear burning. The 
hydrodynamic equations are integrated using the PPM recon 
struction scheme (Colella & Woodward 



questi oned ( |Schneider et al.|1999| |Turkel|1999"||Almgren et al 



2006). However, a recent study by Meakin & Arnett (2007a 
based on a direct comparison of anelastic ( Kuhlen et al.||200? 



solver for real gases according to Colel 



a & Glaz| (T984| . The 



chemical species are evolved by a set of additional coupled con- 



tinuity equations (Plewa & Miiller 1999| ). Source terms in the 
momentum and energy equations due to self-gravity and nuclear 
burning are treated by means of operator splitting at the end 
of every integration step. The gravitational potential is approx- 
imated by a one-dimensional Newtonian potential which is ob- 
tained from the spherically averaged mass distribution. The stiff 
nuclear reaction network is integrated using the semi-implicit 
Bader -Deuflhard method, ( [Bader & Deuflhard|1983)|Press et aT] 
1992 ) which allows for large time integration steps. 

In Herakles, a program cycle consists of two hydrodynamic 
time steps and proceeds as follows (in 2D simulations Step (3) 
is omitted): 

1 . The hydrodynamic equations are integrated in r-direction (r- 
sweep) for one time step including the effects of heat con- 
duction. The time averaged gravitational forces are com- 
puted, and the momentum and the total energy are updated 
to account for the respective source terms. Subsequently, the 
equation of state is called to update the thermodynamic state 
due to the change of the total energy. 

2. Step (1) is repeated in ^-direction (fl-sweep). 

3. Step (1) is repeated in ^-direction (0-sweep). 

4. The nuclear reaction network is solved in all zones with sig- 
nificant nuclear burning (i.e, where T > 10 8 K), and then 
the equation of state is called to update the pressure and the 
temperature. 

5. In the subsequent time step, the order of Step (1), (2), and (3) 
is reversed to guarantee second-order accuracy of the time 
integration, and Step (4) is repeated with the updated state 
quantities. 

6. The size of the time step for the next cycle is determined. 

The velocities in the convection zone, even at the peak of 
the core helium flash, are subsonic corresponding to Mach num- 
bers M ~ 0.01 ( [Mocak et"aL"1|2008| ). In this regime, the appli- 
cability of Riemann solver based methods, like PPM, has been 



1984), and a Riemann [5] 



and fully compressible simulations computed with PPM shows 
that at Mach numbers around 10 -2 , PPM can capture the proper- 
ties of convective flows well. 



4. Hydrodynamic simulations 

We performed two two-dimensional and one three-dimensional 
simulation whose properties are summarized in Table[2] [] All 
simulations, except model hefl.2d.b, cover 6000 s of evolution- 
ary time during the peak of the core helium flash, and were com- 
puted on a computational grid spanning a wedge of 120° in an- 
gular directions centered at the equator. The grid had a radial 
resolution of Ar = 5.55 x 10 6 cm, and an angular resolution of 
AO — A<p = 1.5°). The rather wide angular grid appeared to be 
necessary for the three-dimensional simulations, due to the size 
of the largest vortices (~ 40°) found in previous two-dimensional 
simulations (Moca k et al.|2008| ). 

The 2D model hefl.2d.b was evolved for almost 34 hrs 
(130000 s) on a grid covering the full 180° angle in 9 direc- 
tion. It was simulated at almost twice the resolution as the other 
models. The properties of model hefl.2d.b and its evolution are 
discussed in Section|6j and to some extent also in Sections|4], and 



All models posses a convection zone that spans 1.5 density 
scale heights and that is enclosed by convectively stable layers 
extending out to a radius of 1.2 x 10 9 cm. We used reflective 
boundary conditions in every coordinate direction and for all 
models, except model hefl.2d.a. For this model, it turned out 
to be necessary to impose periodic boundary conditions in an- 
gular direction, because reflective boundaries together with the 
120 ° wedge size affect the large scale convective flow adversely 
leading to higher convection velocities. 

We are limited with our explicit hydrodynamics code by the 
CFL condition which is most restrictive near the center of the 
spherical grid. Therefore, we cut off the inner part of the grid at 
a radius of 2 x 10 8 cm that is still sufficiently far away from the 
radius of the temperature inversion at r ~ 5 x 10 8 cm. To trigger 
convection, we imposed a random flow field with a maximum 
(absolute) velocity of 10 cm s _1 , and a random density pertur- 
bation Ap/p < 10~ 2 . Imposing some explicit non-radial per- 
turbations is necessary, because a spherically symmetric model 
evolved with Herakles on a grid in spherical coordinates will re- 
main spherically symmetric forever. This is different from the 
study of Asida & Arnett ( 2000| ), who did not perturb the initial 
model by an artificial random flow, as instabilities were grow- 
ing from round-off errors. The different perturbation techniques 
seem not to influence the final thermally relaxed steady state 
dMeakin & Arnett|2007b| ). 

Because thermal transport of energy by conduction and ra- 
diation is roughly seven orders of magnitude smaller than the 
convective energy flux, it has been neglected in our simulations. 
Most of the liberated nuclear energy is carried away by con- 
vection. All computed models are non-rotating, because rotation 
seems not to play an important role during the core helium flash 
( |Lattanzio et al.|2006| ). 



We simulated another 3D model at lower resolution than model 
hefl.3d (using a wedge of 90 x 80 x 80 zones centered at the equator 
with Ar = 11.1 x 10 6 cm and A0 = A0 = 1.5°), but we do not discuss 
this model here any further, because of its larger numerical viscosity 
and its almost 50% smaller convective velocities. 
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Table 2. Some properties of the three and two-dimensional simulations: number of grid points in r (N,-), 8 (Ng) and , <p (Afy) direction, 
spatial resolution in r (Ar in 10 6 cm), 8 (A8), and <p (A0) direction, characteristic velocity v c (in 10 6 cm s -1 ) of the flow during the 
first 6000 s, expansion velocity at temperature maximum v exp (in m s ), typical convective turnover time t Q - 2R/v c (in s) where 
R is the height of the convection zone, and maximum evolutionary time t max (in s), respectively. 
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grid 
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A(p 


v c 
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Fig. 2. Temporal evolution of the angular averaged kinetic en- 
ergy density (in erg s ) of models hefl.3d (upper) and hefl.2d.a 
(lower), respectively. 



4. 1 . Temporal evolution 

The 3D model hefl.3d and the corresponding 2D model hefl.2d.a 
evolve initially (t < 1200 s) similarly. Convection sets in af- 
ter roughly 1000 s (the thermal relaxation time), and hot rising 
bubbles appear in the region where helium burns in a thin shell 
(r ~ 5 x 10 8 K). After another ~ 200 s, the bubbles cover the com- 



cn 

10 42 




Time (10 s) 



Fig. 3. Temporal evolution of the total kinetic energy of mod- 
els hefl.3d (solid) and hefl.2d.a (dotted), respectively. 



plete height of the convection zone (R ~ 4.8 x 10 8 cm). The flow 
eventually approaches a quasi-steady state consisting of several 
upstreams (or plumes) of hot gas carrying the liberated nuclear 
energy off the burning region, thereby inhibiting a thermonuclear 
runaway. 

The models exhibit a sandwich-like stratification: two sta- 
ble layers enclose the convection zone at the top and bottom, 
respectively (Fig.[2]l. The convection zone is characterized by a 
large kinetic energy density, and the adjacent convectively stable 
layers show waves induced by convection. The region of high ki- 
netic energy density appearing at ~ 3000 s in the bottom layer in 
model hefl.3d (see top panel of Fig. [2]) is an artifact caused by the 
proximity of the reflective inner radial boundary (at 2 x 10 8 cm) 
of the computational grid. 

The models reach a steady state after roughly 2000 s (Figs. [2] 
and [5J. The steep increase of the total kinetic energy from 
10 39 erg to 10 45 erg (Fig.f) marks the onset of convection. The 
kinetic energy density shows small fluctuations in the fully 
evolved convection zone, and is by an order of magnitude larger 
in model hefl.2d.a than in model hefl.3d. This is in agreement 
with other studies, as it is well know that two-dimensional tur- 
bulence is more intensive (see, e.g., |Muthsam et aL] p995 )). The 
total energy production is about 10% higher in the 3D model due 
to its slightly higher temperature and the strong dependence of 
the triple-a reaction rate on temperature. 

We observe buoyancy driven gravity waves in the convec- 
tively stable layers (Fig.|2"l Fig.[9| (|Zahn||1991| |Hurlburt efaL 



1986, |T994l |Meakin & Ajnett|[2007bl i~but we do not discuss 
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Fig. 4. Upper panels: Snapshots taken at ~ 4754 s showing contour plots of the absolute value of the velocity (in units of 10 6 cm s _1 ) 
for the 2D model hefl.2d.a (left), and for the 3D model hefl.3d in a meridional plane of azimuth angle <p = 50° (middle) and <p = 70° 
(right), respectively. Lower panels: Normalized power spectra of angular fluctuations in the absolute velocity as a function of radius 
for the 2D model helf.2d.a (left) averaged over time from 2100 s to 9500 s, and for the 3D model hefl.3d (right) averaged over time 
from 2250 s to 6000 s, and azimuthal angle. The dashed vertical lines mark the edges of the convective zone of the initial model 
according to the Schwarzschild criterion. 



these waves any further here, because their properties are likely 
biased by the reflective boundaries ( |Asida & Arnett|2000| l. We 
only point out that the differences in amplitude and frequency 
seen in Fig. [2] are physical, i.e., gravity waves have a lower fre- 
quency and amplitude in 3D than in 2D, although the energy 
carried by them is likely similar (Kira ga et al.|1999) . 

4.2. Structure of the convective flow 

Fully evolved convection (t > 2000 s) in the 3D model hefl.3d 
differs significantly from that in the corresponding 2D model 
hefl.2d.a. The convective flow is dominated in the 2D model 
by vortices having (angular) diameters ranging from 30°to 



50°(Fig|4|, and an aspect ratio of close to one. The vortices are 
qualitatively similar to those found in other two-dimensional 
simulations (|Hurlburt et al.||1986| |1994[ |Porter & Woodward| 
|1994[ |Bazan & Arnett||1998| >. This vortex structure of 2D tur- 
bulence is quite typical, and arises from the self-organization of 
the flow ( |Fornberg| 1 977) |Mc Williams 1 1 984) . 

The convective flow in the 3D model hefl.3d consists of 
column-shaped plumes (Fig.|4] and [9]), and contrary to the 2D 
model hefl.2d.a, it does not show any dominant angular mode. 
The typical angular size of turbulent features ranges from 10° to 
30° (Fig|4]i. The power spectra of angular velocity fluctuations 
show that turbulent elements have an almost time-independent 
characteristic angular size of 30° - 50°in case of the 2D model, 
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Table 3. Root mean square fluctuation amplitudes of various 
variables within the convection zone (cnvz: 5 x 10 8 cm < r < 
9.2 x 10 8 cm) averaged over a time period of about 2500 s: tem- 
perature T /(T), density p /(p), helium abundance 4 He /( 4 He), 
and carbon abundance 12 C /( 12 C). 



pos 


run 


T'I(T) 


P/(P) 


4 He'/< 4 He> 


12 c'/< 12 c> 


cnvz 


hefl.3d 
hefl.2d.a 


0.00058 
0.00074 


0.00015 
0.00021 


0.00009 
0.00007 


0.01433 
0.01272 



while the spectra computed for the 3D model change with time 
and exhibit no dominant angular mode. 

We find turbulent flow features across the whole convection 
zone resulting from the interaction of convective up and down- 
flows. Close to the edges of the convective zone we observe the 
smallest turbulent flow features that form when compact turbu- 
lent plumes are decelerated and break-up ( Brummell et al. 2002} . 
Shear instabilities likely play an important role in the develop- 
ment of the turbulent flow as well (Cattaneo et al.|199lj ). 

Tables[3]and 4]provide time-averaged root mean square fluc- 
tuation amplitudes of various variables of the convective flow 
inside the convection zone and near its edges, respectively for 
models hefl.3d and hefl.2d.a. The temperature and density fluc- 
tuations are 30-40 % larger in the 2D model than in the 3D one. 
This is expected, because vortices are stable in 2D flows but de- 
cay in 3D ones (FigQ. The fluctuations in the composition ( 4 He 
and 12 C) are larger by 10-30% in the 3D model which is a re- 
sult of more "broken" and hence more non-uniform mixing of 
chemical elements than in the 2D one. 

The temperature and density fluctuation amplitudes are a fac- 
tor of 2-3 larger near the inner edge than near the outer edge of 
the convective layer. At both edges are the fluctuation amplitudes 
by a factor of 2-4 larger in 2D than in 3D models. The fluctua- 
tions in the composition ( 4 He and 12 C) in both models differ 
close to the convective boundaries as well, by up to a factor of 
two. 



Table 4. Root mean square fluctuation amplitudes of various 
variables at the inner (r — 5 x 10 8 cm) and outer (r = 9.2 x 
10 8 cm) edge of the convection zone averaged over a time pe- 
riod of about 2500 s: temperature T /(T), density p lip), helium 
abundance 4 He /( 4 He), and carbon abundance 12 C /( 12 C). 
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Fig. 5. Auto-correlation function A[(ro;r) measuring the radial 
extent of flow patterns (see Eq.[T]i at different radii ro (4.8 x 
10 s cm, 5.9x 10 8 cm, 7x 10 8 cm, 8.1 x 10 8 cm, and 9.25x 10 8 cm) 
and t ~ 4000 s measuring the radial extent of flow patterns for 
the 2D model hefl.2d.a (dotted lines) and the 3D hefl.3d (solid 
lines), respectively. 
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run 


T'KT) 


plip) 


4 He'/( 4 He> 


12 c'/< 12 c> 




hefl.3d 


0.00643 
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Fig. 6. Auto-correlation function A2(ro; t—to) measuring the life- 
time of flow patterns (see Eq.[2]) at two different epochs for the 
3D model hefl.3d (solid: f ~ 2260 s, dashed: t ~ 2900 s), and 
for the 2D model hefl.2d.a (crosses, diamonds), respectively. The 
radius ro is 7.6 X 10 8 cm. 



4.3. Stability of flow structures 

To analyze the size and the stability of the vortices we introduce 
an auto-correlation function of the radial velocity 



Ai(r ;r) = 



(v r (r ) v r (r))n, t 



<v,.(ro) 2 )^<v V (r) 2 )^; 



2xV2 



(1) 



that measures the radial extent of flow patterns at a given ra- 
dius ro, or equivalently the radial size of vortices. The notation 



()n,r indicates averaging over angles and time. A second auto- 
correlation function 



A 2 (r ; t - t ) = 



(v r (ro, to) v r (rp, t)) n 
<v r (ro,/o) 2 >o /2 <v,(r ,f) 2 )o /2 



(2) 



provides a measure of the lifetime of flow patterns at radius 
ro, beyond time f . Here, we average only over angles. Both 
correlation functions have the properties -1 < A\ ^ < 1, and 
Ai(ro\ ro) = 1 and A 2 (ro; to - to) = 1, respectively. They are sim- 
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ilar to the autocorrelation function used by |Chan et al.| ( [l982) ; 
|Chan & Sofia1 ( |T986) . 

Figure |5J displays the radial auto-correlation for models 
hefi.2d.a and hefl.3d.a and confirms the extension of the con- 
vective flow across the whole convective region as determined 
by the Schwarzschild criterion in the initial stellar model. The 
broad plateaus with A\ as 1 corresponding to the 2D model 
hefl.2d.a bear eveidence of the axial symetry imposed in the two- 
dimensional case which leads to pronouced circular vortices. In 
the three-dimensional case, the distributions of A\ tend to differ 
from unity at nearly all radii. 

Figure[6] shows the temporal auto-correlations with two typ- 
ical results for different to- The three-dimensional model always 
shows the typical behaviour of a decrease of the function value to 
0.5 within 200-250 s. From this we conclude that the flow pattern 
fluctuates always in the same way. This is different in the two- 
dimensional model, where A2 can keep high values for a much 
longer time implying rather persistent structure (the vortices) of 
the convective flow. 



6 




4.4. Turbulent entrainment and the width of the convection 
zone 

Convection may induce mixing in convectively stable layers 
adjacent to convectively unstable regions. Following Meakin 
|& Arnett| ( |2007b| l, we prefer to call this process turbulent en- 
trainment (or mixing), a term also known in oceanography, see 
e.g., Fernando (1991). The commonly used term (convective) 
overshooting accounts only for localized ascending or descend- 
ing plumes crossing the edge of the convection zone. If the fill- 
ing factor of these plumes or their crossing frequency is high, 
they can change the entropy in convectively stable regions sur- 
rounding convection zones, a process that is known as penetra- 
tion (Bmmmell et al. 2002). Turbulent entrainment accounts for 
both overshooting and pe netra t ion. 

Contrary to Hurlburt et al. ( 1994 1 we determine the depth of 
the entrainment neither by the radius where the kinetic energy 
carried into the stable layers is zero, nor by the radius where the 
kinetic energy has dropped to a certain fraction of its maximum 
value (Brum mell et al.| |2002). We find both conditions insuffi- 
cient, because the kinetic energy flux becomes zero much faster 
than the convective flux, which is another possible indicator of 
the depth of the entrainment (see next subsection). Instead, we 
rather use the 12 C mass fraction, as it is low outside the convec- 
tion zone during the flash (X( 12 C) < 2 x 10~ 3 ), and as it can rise 
there only due to turbulent entrainment. 

In this study, we use the condition X( 12 C) = 2 x 10~ 3 to 
define the edges of the convection zone. Due to the turbulent 
entrainment, these edges are pushed towards the stellar center 
and surface (Fig.|7]i. Hence, the width of the convection zone in- 
creases on the dynamic timescale, which is in contradiction with 
the predictions of one-dimensional hydrostatic stellar modeling, 
where the width of the convection zone is determined by the lo- 
cal Schwarzschild or Ledoux criterion. However, the criterion 
for the width of the convection zone cannot be a local one due to 
turbulent entrainment caused by convection. 

The speed, at which the radius of the outer edge of the con- 
vection zone increases with time due to entrainment, is estimated 
for models hefl.2d.a and hefl.3d to be at most 14 m s . The ra- 
dius of the inner edge of the convection zone changes at a much 
smaller rate ( [Bazan & Arnett|1998||Meakin & Arnett|2007b) l, as 
the region interior to the convection zone is more stable against 
convection and has a higher density than the region exterior to 
the convection zone ( |Singh et al.|1 995). The entrainment at the 
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Fig. 7. Carbon mass fraction X( 12 C) (top), entropy S (middle), 
and entropy gradient VS (bottom) as a function of radius near 
the outer edge of the convection zone of model hefl.3d at three 
different epochs: 1 1 = 2000 s (dashed), 1 1 = 4000 s (dash-dotted), 
and fi = 6000 (solid). In addition, the initial X( 12 C) profile is 
shown in the top panel (dotted). 
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Fig. 8. Square of the Brunt-Vaisala buoyancy frequency as a 
function of radius for the 3D model hefi.3d at t = 4638 s. 



bottom of the convection zone also leads to a heating of the cool 
interior layers ( jDeupree & Coie][T983j |Cole et al.|[T985| ). This 
seems to be a robust feature of convection zones driven by nu- 
clear burning, and is observed in other studies too (e.g., Asida & 
ArnettlpOOOl l). 

The region just interior to the convection zone shows less 
entrainment ( |Bazan & Arnett|[T998l |Meakin & Amett||2007b] l, 
as the square of the Brunt- Vaisala frequency is almost ten times 
larger there than that in the region just outside the outer edge 
of the convection zone. The Brunt-Vaisala frequency is a good 
stability indicator since it is related to the behavior of convec- 
tive elements within a convection zone, a fact also pointed out 
by Hurlburt et al. ( 1994 1. The Brunt- Vaisala frequency can be 
written as ( TMeakin & Arnett|2007b| l : 



N 



d\np <91np 

~d~r 



d) 



(3) 



where g, p ,and r are the gravitational acceleration, the den- 
sity and the radius, respectively. A layer is con vectively stable, 
if N 2 > 0, and unstable otherwise [^] (see, e.g., Kippenhahn & 
Weigert| ( [T990| ). 

The Brunt-Vaisala frequency differs in the 2D and 3D simu- 
lations slightly, and has on average a very small negative value 
inside the convection zone: N 2 = -4.4 x 10 6 rad 2 crrT 2 in the 
2D model hefl.2d.a, and N 2 = -1.2 x 1(T 6 rad 2 cnr 2 in the 3D 
model hefi.3d. Just outside (inside) the inner (outer) edge of the 
convection zone of the 2D model hefi.2d.a we find N 2 = 0.580 
rad 2 cm" 2 (0.053 rad 2 cnr 2 ), andN 2 = 0.583 rad 2 cnr 2 (0.052 
rad 2 crrT 2 ) for the 3D model hefi.3d, respectively (Fig. 8]). In 
the 2D model hefl.2d.b, these frequencies are higher by about a 
factor of two. 

Due to its high stability, the radius of the bottom edge of 
the convection zone did not change during the time covered by 
our simulations, except for an initial jump over one radial grid 
zone from r = 4.69 x 10 8 cm to r = 4.63 x 10 8 cm, when it was 



2 The Brunt-Vaisala frequency is related to the bulk Richardson num- 
ber known from oceanography, which is a more direct measure of the 
stability of the edges of a convection zone in the presence of a turbulent 
flow dMeakin & Arnett|2007bli. 



touched by convective downflows for the first time. However, en- 
trainment may move the edge further towards the stellar center 
later in the evolution (see Sect.|6]l. The entrainment rate (i.e., the 
velocity at which the convective boundary moves) is lower by a 
factor of ~ 5-6 at the bottom of the convection zone in our 2D 
simulations of the core helium flash as compared to that at the 
outer boundary (~ 14 m s ; see above). This behavior was also 



observed in 3D simulations of oxygen burning shell (Meakin & 



Arnett 2007m). This implies that the entrainment rate at the inner 
edge of the convection zone is ~ 2.5 m s in our core helium 
flash simulations. The corresponding change in radius is only 
1 .5 x 10 6 cm or about a quarter of the width of a radial zone dur- 
ing the time covered by the simulations, and hence too small to 
be seen. As these estimates are resolution dependent, the values 
presented should be considered as order of magnitude estimates, 
only. 

The entrainment is more efficient in the 2D model hefl.2d.a 
than in the 3D model hefl.3da. In 2D, the observed convective 
flow structures are large and fast rotating vortices that due to the 
imposed axisymmetry are actually tori (B azan & Arnett|1998| l. 
They have a high filling factor near the edge of the convec- 



tion zone where they overshoot or penetrate (Brummell et al. 
2002 ). 3D structures crossing the edge of the convection zones 
are smaller (localized) plumes with a lower filling factor and 
smaller velocities. 

We studied convective stability in detail only for the layer 
above the convection zone, where the gas is only weakly degen- 
erate, and we are thus able to compare our results with those 
of similar systems simulated in 3D at high Reynolds numbers 
(~ 10 4 ) by Brummell et al. (2002 1. According to these authors, 
a stable layer allows for more "overshooting", if its entropy gra- 
dient is smaller. We found that the turbulent entrainment lowers 
the entropy gradient at the outer edge of the convection zone in 
our simulations (Fig. [7]). Hence, the stability of the exterior stable 
layer decreases with time, allowing for more turbulent entrain- 
ment. 

Contrary to Brummell et al. ( 2002[ l, who studied only single 
fluid flow, our simulations involve a mixture of different fluids 
of different composition. This complicates the above argumen- 
tation, as in a multi-fluid flow a shallow entropy gradient does 
not necessarily imply that the layer is less stable against turbu- 
lent entrainment. We plan to address this issue in more detail 
elsewhere. We have not analyzed the stability of the region be- 
low the inner edge of the convection zone, because it is highly 
degenerate and appears to be very stable. No significant entrain- 
ment of 12 C was observed there during the entire simulation of 
models hefl.2d.a and hefl.3d, respectively. Simulations covering 
longer periods of time are therefore needed (see Sect. [6]). The 
analysis of the energy fluxes inside the convection zone and near 
the outer edge of the convection zone, which will be discussed 
in the next subsection, will provide more information about the 
phenomenon of the turbulent entrainment. 



4.5. Energy fluxes 

Energy fluxes are a useful tool for understanding convective 
flows ([Chan & Sofia|1986l|Hurlburt et al.|1986||Muthsam et aL 
[19951 Muriburt et al.||1994j |Brummell et al.||2002| |Meakin & 
|Arnett|2007b I. They allow one to discriminate the energy trans- 
port due to different processes and mechanisms (e.g., due to con- 
vection, heat conduction, etc). We thus analyzed various energy 
fluxes and source terms that are defined in Appendix[B| All en- 
ergy fluxes and source terms (see, Fig.[T0]> are averaged over sev- 
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Fig. 9. Different views of isosurfaces of the velocity field for the 3D model hefl.3d at t — 3000 s. The blue isosurface corresponds to 
radial downflows with v r — -6 x 10 5 cm s _1 , and the yellow and brown isosurfaces show radial upflows with v r — 6 x 10 5 cm s _1 , 
and 1 x 10 4 cm s _1 (gravity waves), respectively. The edge sizes of the box are 1.2 x 10 9 cm and 2.4 x 10 9 cm, respectively. The 
yellow-greenish sphere in the bottom right panel marks the top of the convection zone according to the Schwarzschild criterion. 



eral convective turnover times, as they change considerably with 
time due to the appearance of plumes (Muth sam et al.[ 1995). 

The convective energy flux is mainly positive as heat is 
mostly transported upwards by convection. It is larger in the 
2D model hefl.2d.a, as in 2D the convective flow structures are 
more laminar and ordered, and thus experience less dissipation. 
The smaller value of Fc in the 3D model hefl.3d is a result of a 
less ordered flow throughout the whole convection zone. This is 
in agreement with Table|3] which shows that the fluctuations in 
temperature and density are smaller in the 3D model. 

A kinetic energy flux arises from deviations from the mean 
convective flow {i.e., mainly from the upflow-downflow asym- 
metry). Typically it is largest in the most turbulent regions of 
the flow, where on top of the kinetic flux due to the up- and 
down-flow asymmetry, there is also a significant contribution 



due to the asymmetry of localized turbulent (convective) ele- 
ments. Directly connected to this is the offset between the max- 
ima of the kinetic and the convective flux (see Fig. 10 1, which 
reflects the fact that the convective flow decays more efficiently 
when the flow becomes strongly turbulent. 

The work done by buoyancy forces is positive in the whole 
region of dominant nuclear burning (see Fig. 10 1 indicating ei- 
ther less dense and hot gas moving upwards or more dense and 
cooler gas moving downwards. A negative value of the work 
done by buoyancy forces implies the opposite situation, i.e., less 
dense and hot gas flowing downward or denser and cooler gas 
flowing upward. The latter is known as buoyancy breaking lead- 
ing to a deceleration of the flow motion (Brumm ell et al.|2 002), 
and to the unusual situation that hot matter tends to sink, and 
cool matter is likely to rise. 
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Fig. 10. Various energy fluxes and source terms as a function of radius averaged (from 2000 s to 6000 s) over about four convective 
turnover times. Panel (a) shows the convective flux Fc of the 2D model hefl.2d.a (solid-red) and the 3D model hefl.3d (dashed- 
black) together with the kinetic energy flux F% in the 2D (dash-dotted-red) and 3D (dash-dot-dotted-black) model, respectively. The 
dotted vertical lines mark the edges of the convection zone in the initial model according to the Schwarzschild criterion. Panel (b) 
gives the source terms due to the work done by buoyancy forces Pa (dash-dot-dotted black) and due to volume changes Pp (dashed 
black) in the 3D model hefl.3d, and in the 2D model hefl.2d.a (dash-dotted-red, solid red) respectively. The vertical lines enclose 
the nuclear burning zone (T > 10 8 K). Panels (c) and (d) show an enlarged view of the energy fluxes and source terms displayed in 
panels (a) and (b) near the outer edge of the convection zone. 



The gas should expand while rising up through the con- 
vection zone, i.e., the work done by buoyancy {Pa) and by 
volume changes (Pp) should always be anti-correlated. The 
anti-correlation is clearly seen only in the 2D model hefl.2d.a, 
whereas in the 3D model hefl.3d the quantities are on average 
anti-correlated only in the central region of the convection zone. 
At the inner edge of the convection zone, buoyancy drives gas 
upwards which is on average simultaneously compressed, prob- 
ably by the broad downflows. At the outer edge where buoyancy 
braking occurs, the gas on average expands. Hence, it must ex- 
pand faster than the upflows cool and are being compressed (a 
situation observed for the 2D model hefl.2d.a). 



All the fluxes discussed here agree qualitatively well with 



those of our previous high-resolution 2D simulations (Mocak 
|et al.|2008] >. They are also qualitatively very similar to those of 
the high-resolution 3D simulations of Brummell et al. (2002) 
who investigated a stratified model with a convectively stable 
region located on top of a convectively unstable region both 
consisting of an ideal gas with a very high Reynolds number 
(~ 10 4 ). The angular and time averaged radial distributions of 
the kinetic and convective fluxes seem to be robust in the con- 
vectively unstable region, as our 3D results are qualitatively sim- 
ilar to those obtained in sev eral other 3D studies (|Hurlburt et"aL] 
[T994l|Brummell et al.|2002"l|Meakin & Arnett|2007b] i~ 
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The outer part of the convection zone, where buoyancy 
breaking occurs, resembles the overshooting region due to ac- 



tive penetration of plumes described by |Hurlburt et al. (1986 
1 1 994> and |Brummell et al.| ( |2002| >. Note, however, that this re- 
gion is convectively stable at the beginning of their simulations, 
i.e., buoyancy breaking takes already place inside the convection 
zone in our models. 



2.0 



The distribution of the kinetic flux (Fig. 10 1 exhibits struc- 
tural differences in the convection zone between 2D and 3D 
flows. The typical evolved 2D flow contains well defined vor- 
tices (Fig. [4]) whose central parts never interact with the region 
of dominant nuclear burning. This results in a reduced kinetic 
energy flux at r ~ 5.5 x 10 8 cm, as this region corresponds to 
the central region of the vortices, which do not experience any 
strong radial motion. On the other hand, the distribution of the 
kinetic energy flux in the 3D model hefl.3da is rather smooth as 
a result of the column-shaped flow structures (Fig. [9]). 

The convective flux changes sign in stable layers since the 
downflows or upflows penetrating into the stable zones are sud- 
denly too hot or cold compared to the surrounding gas. The pene- 
tration continues until the momentum of the convective elements 
is used up or diffusion smooths out the perturbations, and the 
convective flux approaches zero (Brummell et al. 2002). This 
fingerprint of turbulent entrainment is clearly present at both 
convective boundaries of our models (Fig. 10 1. The convective 
energy flux is relatively strong even in regions where the kinetic 
energy flux is almost zero. Therefore, the "zero" kinetic flux cri- 
terion (Hurlbur t et al.|199 4; Brummell et al. 2002) seems to be a 
bad indicator of entrainment which may extend well beyond the 
location where the kinetic flux becomes small. 

In fact, what is happening at the convective boundaries is an 
exchange between the potential energy of the stratification given 
by the buoyancy jump db = A^ 2 dr and the kinetic energy of the 
turbulence. Turbulence looses its kinetic energy by doing work 
against gravity, which leads to a reduction of the buoyancy jump, 
and hence to stability weakening of the stable boundary layer to 
the effects of the turbulent entrainment. The buoyancy jump db 
is a direct measure of the stability of the boundary layer. To mix 
gas into the boundary layer the buoyancy must be reduced. This 
is accomplished through the buoyancy flux q = Pa/p, which 
is related to the temporal variation of the buoyancy jump by 
db/dt = -div(q), where Pa and p are the sink/source term of 
the kinetic energy due to buoyancy forces and the density, re- 
spectively. 

The convective flux can directly be related to the buoyancy 
flux that is a function of Pa by a linear relation Fq = F~c(q) de- 
scribed in more detail by Meakin & Arnetf (2007b]). Figure 10 
shows that this agrees well with what we observe in our simula- 
tions. It also supports our previous conclusion that entrainment 
is well indicated by the convective flux and the related buoy- 
ancy flux which via the equation db/dt = -div(q) leads to the 
decrease of the buoyancy jump in the stable layer. This in turn 
reduces the convective stability of that layer. 

The properties of the entrainment at the outer convective 
boundary differs in models hefl.2d.a and hefl.3d (see Fig.[T0|i. 
In the 2D model entrainment is reaching deeper into the stable 
layer due to a more active convection zone with higher typical 
velocities (Fig. 1 1 1 than in the 3D model. The radial distribu- 
tion of the work done by buoyancy Pa is qualitatively similar in 
both the 2D and 3D models, i.e., it is negative indicating buoy- 
ancy breaking. Nevertheless, the work done by gas compression 
or expansion Pp is different. In the 2D model hefl.2d.a the gas 
on average is compressed (Pp is positive), while in the 3D model 
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Fig. 12. Radial distributions of the time (from 2800 s to 4800 s) 
and angle-averaged density (upper) and temperature (lower) 
fluctuations in the 3D model hefl.3d (solid red line). The pan- 
els also show the angular variation of the respective quantity at 
a given radius (gray shaded region). The dotted vertical lines 
mark the edges of the convection zone as determined by the 
Schwarzschild criterion. 



hefl.3d the gas is expanding at the boundary. This again confirms 
that 2D and 3D convective flows are qualitatively different. 

4.6. The flow within the convection zone 

The amount of energy (F c + F K ) which has to be transported 
by convection in order to prevent a thermonuclear runaway dur- 
ing the flash is similar in models hefl.2d.a and hefl.3d. Since the 
convective flux is almost the same in both models, the resulting 
typical convective velocities are higher in the 2D model than in 
the 3D one (Fig. [11} . 

The velocities in the 3D model hefl.3d match those predicted 
by mixing-length theory better than in the 2D model hefl.2d.a, 
where they are about a factor of 2 larger. This behavior was 
also observed in other hydrodynamic simulations of convective 
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Fig. 11. Radial distributions of the time (from 2000 s to 4000 s) and angle-averaged velocity components (v r , vg, v^) and velocity 
modulus (Vabs) for the 2D model hefl.2d.a (left), and the 3D model hefi.3d (right), respectively. The panels also show the velocity 
predicted by the mixing-length theory (v m / f ) . 
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Fig. 13. Fractional volume occupied by upflow and downfiow 
streams averaged over ~2000 s as a function of radius for the 
3D model hefi.3d (solid), the 2D low-resolution model hefi.2d.a 
(dashed-dotted), and the 2D high-resolution model hefl.2d.b 
(dotted), respectively. 



flows; see e.g., |Muthsam et al.| ( [l99^ ; |Meakin & Arnett| ( |2007b| l. 
The radial velocities in the regions above and below the convec- 
tion zone are smaller than the angular velocities in both mod- 
els, which is a typical feature of gravity waves (Asida & Arnett 
[2DD0J. 

In the 3D model hefl.3d, the flow in the convectively sta- 
ble layer beneath the convection zone exhibits some numerical 
artefact's due to the proximity of the inner grid boundary (see 
Fig. Hi. The radial distributions of the time and angle-averaged 



a sharp cut-off at r — 2.5 x 10 8 cm (Fig. 1 1 Fig. 12 1. The sharp 
cut-off is caused by the artificial damping we had to apply to the 
velocity field in the innermost grid zones to prevent numerical 
instabilities from spreading limitless to the convection zone. 

Although, the flow velocities in the 3D model match those 
predicted by mixing-length theory very well, one should keep 
in mind that with increasing resolution the flow velocities will 
likely increase due to the reduced numerical viscosity. This trend 
is confirmed by the velocities obtained for the high-resolution 
2D model hefl.2d.b that are a factor of two higher than in the 
low-resolution model hefl.2d.a. 

Near both edges of the convective zone there are large nar- 
row peaks visible in the radial distributions of the time and 
angle-averaged density and temperature fluctuations (Fig. 12 1. 



These peaks are not caused by compression or expansion, but 
they are a result of the density and temperature discontinuities 
at the edges of the convection zone ( [Meakin & Arne tt 2007a; 
Arnett et al. 2007 ), because any angle-dependent radial pertur- 



components of the velocity field and of the density and tempera- 
ture fluctuations show pronounced maxima at ~ 3 x 10 8 cm and 



bation will cause large angular variations of density and temper- 
ature at these discontinuities. 

The temperature fluctuations within the convection zone are 
rather uniformly distributed, but they are more intense near the 
outer edge of the convection zone, where they are only weakly 
correlated with the radial velocity (Fig.[l4|. At the top of the 
convection zone the emerging rising plumes are embedded in 
an environment which sinks down (Fig. 14 bottom panels). This 
situation is similar to the sinking down-drafts with upwelling 
centers found in simulations by |Nordlund & Dravins ( 199Q] > and 
ICattaneo et aTlfiWfl l. 



4. 7. Upflow-downflow asymmetry 

The 2D and 3D simulations share the common property of an 
upflow-downflow asymmetry (Hurlburt et al. 1994| Muthsam 
|et al.|1993)|Brummell et al.|2002) >. The downflows cover a much 
larger volume fraction of the convection zone than the upflows 
(Fig. 13 1. The filling factor of the downflows increases with de- 
creasing depth ( |Rast et al.|1993||Meakin & Arnett|2007b| l across 
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most of the convection zone, and the downfiows are more dom- 
inant in the 3D model hefl.3d than in the 2D models. Contrary 
to the simulations of oxygen burning shell of Asida & Arnett 



(2000) we find that the absolute velocities are about 40% higher 
in the upflows than in the downfiows. Hence, the downfiows in 
the convection zone of a star at the peak of the core helium flash 
are slower and broader than the faster and narrower upflows. 



4.8. Mixing 

Cuts through the 3D model hefl.3d at t = 4815 s showing the 
angular variation of temperature, 12 C mass fraction, and radial 
velocity at three different radii (Fig.[T4|) demonstrate that the he- 
lium core at the peak of the core helium flash is a very turbulent 
environment at all heights of the convection zone. 

The bottom of the convection zone contains hot filaments 
of gas where the temperature exceeds that of the environment 
by about 1%. The filaments contain ashes from helium burning, 
i.e., I2 C and 16 0, and they move across the whole bottom of the 
convection zone in a random way. The filaments are correlated 
with upflows, as the hot gas of burned matter is forced by buoy- 
ancy to rise towards the top of the convection zone. 

The apparent turbulent nature of the convective flow indi- 
cated by our simulations implies that the treatment of mixing in 
stars as a diffusive process may lead to inaccurate or even incor- 
rect results. Convective flows are rather advective, as suggested 
by |W(X)dwardetaLl ( |2008) . 



5. Mixing length theory and simulations 

Mixing length theory (or MLT) commonly used for treating con- 
vection in stellar evolutionary calculations relies on assumptions 
and parameters that are often chosen based on convenient ad-hoc 
arguments about the convective flow, like e.g., the value of the 
mixing length, the amount of upflow-downflow symmetry or the 
position where, within the convection zone, convective elements 
start to rise ( |Kippenhahn & Weigert|1990| [Weiss et al.|20Q4| >. 

MLT assumes that the temperature of a convective element 
(blob) is the same as that of the ambient medium surrounding 
it when it starts to rise. However, as a blob will not rise until it 
is hotter than the surroundings, this MLT assumption is contra- 
dictory. MLT further assumes that once the blobs begin to rise 
they carry their surplus of heat lossless over a distance given by 
the mixing length before they release it to the surrounding gas 
instantaneously at the end of their path. These assumptions are 
also not fulfilled in general. Our simulations show that convec- 
tive elements typically start their rise deep inside the star from 
the region of dominant nuclear burning where they are acceler- 
ated by buoyant forces. The assumptions of MLT that convective 
blobs form and begin their motion at different depths of the con- 
vection zone, and that the average convective blob propagates a 
distance equal to half of the assumed mixing length before dis- 
solving with the surrounding gas (Kippenhahn & Weigert 1990), 
therefore do not hold. MLT finally also assumes a correlation be- 
tween the thermodynamic variables and the velocity of the flow 
in a convection zone. However, the results of our simulations 
falsify this assumption (Fig.[l4j). 

According to MLT, the temperature fluctuations in a con- 
vection zone are directly proportional to the mixing length 
and to the deviation of the temperature gradient of the model 
V sml = (d In T Id In p) s [ m from the adiabatic one V a( j = 
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Fig. 15. Upper panel: Radial distributions of the adiabatic tem- 
perature gradient V a d (dotted) and of the temperature gradients 
of model hefl.3d (solid), respectively. The latter distribution is 
a linear fit to the gradients averaged over angle and over the 
first 200 s of the evolution of the model. The gray shaded region 
marks the convection zone CVZ. Lower panel: Same as above, 
but showing the radial distributions in the evolved convection 
zone averaged over roughly 3000 s of evolutionary time. 



(d \nT/d lnp) a( j: 



T 

Y 



^sim ^ad-* 



1 A 

~H p 2 



(4) 



where A is the mixing length, H p the pressure scale height, T 
the absolute value of the temperature deviation from the mean 
(horizontally averaged) temperature T, and p the pressure, re- 
spectively. 

Since T' /T and V sml - V a( j can directly be obtained from 
our simulations, we attempted to test MLT in a qualitative man- 
ner. Our simulations show that in the outer part of the convec- 
tion zone, i.e., in the region where the buoyancy force is getting 
smaller (Fig.[l0| panel b), the temperature gradient of the mod- 
els Voim becomes sub-adiabatic (see lower panel of Fig. 15 I. 



Equation Q which was derived for the adiabatic rise of convec 
tive bubbles would then imply that the temperature of convec 
tive elements should be lower than that of the surrounding gas 
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Fig. 14. Cuts through the 3D model hefl.3d at t = 4815 s showing the angular variation of temperature (in units of 10 s K; left 
panels), 12 C mass fraction (in units of 10~ 3 ; middle panels), and radial velocity (in units of 10 s cm s _1 ; right panels), respectively, 
at three different radii: r\ = 4.8 x 10 8 cm (temperature maximum; top), r-i - 6.5 x 10 8 cm (center of the convection zone; middle), 
and - 9.3 x 10 8 cm (top of convection zone, bottom). The black lines in the right panels mark the boundaries between positive 
and negative radial velocities. 



(hence no convection) in the outer part of the convection zone, or 
that the value of A should be negative. However, convective ele- 
ments do not rise adiabatically in our hydrodynamic simulations 
and the sub-adiabatic gradient means only that the convective 
elements start to cool faster than their surroundings. It does not 
imply necessarily that the elements are already cooler than the 
surrounding gas which would prevent the gas from being con- 
vectively active. Note that initially, the temperature gradient is 
super-adiabatic in the whole convection zone (see upper panel 
of Fig.[T5|, because the stellar evolutionary model used as initial 
input for our simulations is computed under the assumptions of 
the MLT. 



6. Long-term evolution 

Two-dimensional simulations are biased due to the imposed 
symmetry restriction which tend to overestimate the activity in 
the convection zone, but qualitative similarities with geometri- 
cally unconstrained 3D models exists. Moreover, many phenom- 
ena that we observe in 2D models happen also in the 3D ones, 
just slower. Hence, we think it is justified to explore the long- 
term evolution (i.e., covering a few tens of hours instead of a 
few hours) of our initial core helium flash models by performing 
cost-effective 2D instead of very costly 3D simulations. 

In the following we describe the long-term evolution of a 
2D model whose early evolution, covering 8 hrs, was discussed 
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Fig. 16. Snapshots of the spatial distribution of the velocity modulus |v| (in units of 10 6 cm s 1 ) for the 2D model hefi.2d.b at 24 000 s 
(left), 60000 s (middle), and 120000 s (right), respectively. 



in detail by Mocak et al. ( 2008[ ). The model is characterized by 
a very dynamic flow involving typical convective velocities of 
1.8 x 10 6 cms -1 . Our long-term hydrodynamic simulation of 
this model covering 36 firs (see Fig.[T7j) has revealed that the 
global and angle-averaged maximum temperatures continue to 
rise at the initial rate of 40 K s" 1 that is 60% lower than the 
rate predicted by stellar evolutionary calculations. As a conse- 
quence, the typical convective velocities increase by about 50% 
and reach a level of 2.8 x 10 6 cm s _1 at the end of the simulation 



(Fig.[T§. 

Hydrodynamic simulations of convection driven by nuclear 
burning covering several convective turnover times show a rapid 
growth of the convection zone due to the turbulent entrainment 
(Meakin & Arnett 2007b ). An analysis of our simulations based 
on a tracing of the radial position of the convective bou ndaries 
(defined by the condition X( 12 C) = 2xl0~ 3 ; see Sect.|45jl, shows 
a similar behavior. 

Turbulent motion near the upper edge of the convection zone 
pumps material into the convectively stable layer at an entrain- 
ment rate of 14 m s without any significant slowdown over the 
whole duration of the simulation (Fig. 18 1 covering ~ 130000 s 
(or more than 250 convective turnover times). The entrainment 
rate at the inner convective boundary is about a factor of six 
smaller (2.3 m s _1 ) slightly increasing during the second half of 



the simulation (f > 60000 ; see Fig. 18 1. These entrainment rates 



have to be considered as upper limits because of the imposed 
axisymmetry which leads to exaggerated convective velocities 
and large filling factors for the penetrating plumes. The turbu- 
lent entrainment causes a growth of the convection zone on a 
dynamic timescale, in agreement with the oxygen shell burning 
hydrodynamic models of Meakin & Arnett ( 2007b| l. 

Both interfaces at the edges of the convection zone remain 
sharp during the whole length of the simulation (Fig.[T8]l. The 
entrainment is correlated with a decrease of the temperature at 
the outer edge of the convection zone due to the decreasing en- 
tropy (Fig.|7]i. At the inner edge of the convection zone the en- 



00 

O 




2 4 6 8 10 12 14 
Time ( 1 4 s) 

Fig. 17. Temporal evolution of the horizontally averaged tem- 
perature maximum {T) max (solid), and of the global temperature 
maximum T max (solid thin) in the long-term 2D model hefl.2d.b. 
The dashed line correspond to the temporal evolution of the max- 
imum temperature in the stellar evolutionary calculation. 



trainment leads to heating of the cold region at r < 5 x 10 8 cm 
that is cooled by neutrinos (Fig. 18 i. Contrary to the finding of 



|Asida & Arnett| ( |2000[ l, the heating does not penetrate deeper 
into the star than the mixing of 12 C and the other nuclear ashes. 

Due to the growth of the convection zone and due to nuclear 
burning the mean 12 C, 4 He, and ie O mass fractions change in 
the convection zone at rates listed in Table|5] The mean value 
of the I2 C mass fraction in the convection zone decreases at a 
rate of -7.4 x lO" 9 *" 1 until t = 40000 to a value of X( 12 C) 
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Fig. 18. Upper panels: Angular averaged 12 C mass fraction (dashed) and temperature stratification (thick) as a function of radius 
near the inner (left) and outer edge (right) of the convection zone in the long-term 2D model hefi.2d.b at t — 60000 s and t = 
120000 s, respectively. The vertical dotted lines mark the boundaries of the convection zone at t = s. Lower panels: Temporal 
evolution of the position of the inner (left) and outer (right) edge of the convection zone in model hefl.2d.b, respectively. 



Table 5. Approximate rates at which the mean mass fractions of 
4 He, I2 C, and 16 evolve in the long-term 2D model hefi.2d.b 
in the convection zone within the first 40000 s (/?,; in units of 
10~ 9 ), and within the time interval 40000 s to 130000 s (R f ; in 
units of 10~ 9 ), respectively. The quantities X t and Xf give the 
initial (f = s) and final (t = 130 000 s) mass fraction of 4 He and 
mass fractions of 12 C, and ie O abundances (in units of 10" 3 ), 
respectively. 



element 




R f 


X, 


*/ 


4 He 


+7.74 


-7.85 


0.975 


0.974 


12 C 


-7.41 


+7.71 


5.502 


5.918 


16Q 


-0.58 


+ 1.13 


0.927 


1.008 



= 5.2 x 10 3 . Then it begins to increase again at roughly the 



same (absolute) rate +7.7 x 10~ 9 s~' . The 12 C mass fraction de- 
creases, because the volume of the convection zone grows ini- 
tially almost discontinuously due to the sudden start of the en- 
trainment. Hence, nuclear reactions are for a start unable to pro- 
duce enough carbon to compensate for the volume increase. At 
t ~ 110000 s the 12 C mass fraction has risen again to its initial 
value of X( 12 C)= 5.5 X 10~ 3 , and at the end of the simulation at 
t = 130000 s the I2 C mass fraction is 5.9 x 10~ 3 , a value that is 
7% higher than the initial one. 



The mean le O mass fraction shows a similar trend as that of 
12 C as its production depends directly on the 12 C mass fraction. 
The mass fraction of 4 He rises within the first 40 000 s because 
convection is dredging up fresh 4 He from the convectively stable 
layers. Later in the evolution the mass fraction of 4 He decreases, 
as it is being constantly burned. 
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7. Summary 

We performed, analyzed and compared 2D axisymmetric and 3D 
hydrodynamic simulations of the core helium flash. In agree- 
ment with our previous study of 2D hydrodynamic models of 
the core helium flash we find that the core helium flash in three 
dimensions neither rips the star apart, nor that it significantly al- 
ters its structure. 

The evolved convection of the 3D models differs qualita- 
tively from that of the axisymmetric ones. The typical convec- 
tive structure in the 2D simulations is a vortex with a diameter 
roughly equal to the width of convection zone, whereas the 3D 
structures are smaller in extent and have a plume-like shape. The 
typical convective velocities are much higher in the 2D mod- 
els than in the 3D ones. In the latter models the convective ve- 
locities tend to fit those predicted by the mixing length theory 
better. Both 2D and 3D models are characterized by an upflow- 
downflow asymmetry, where downflows dominate. 

Our hydrodynamic simulations show the presence of turbu- 
lent entrainment at both the inner and outer edge of the convec- 
tion zone, which results in a growth of the convection zone on a 
dynamic timescale. While entrainment occurs at an almost con- 
stant speed at the outer boundary of the convection zone, it tends 
to accelerate at the inner one. The entrainment rates are higher 
at the outer edge than at the inner edge of the convection zone, 
as the latter is more stable against entrainment. 

The upper part of the evolved convection zone is character- 
ized by a sub-adiabatic temperature gradient, where buoyancy 
breaking takes places, i.e., rising convective plumes start to slow 
down in this region and eventually descend back into deeper lay- 
ers of the star. Convection should not exist in that mass layer 
according to mixing length theory. 

The fast growth of the convection zone due to entrainment 
has some potentially interesting implications. As entrainment 
is not considered in canonical stellar evolutionary calculations, 
stars evolving towards the core helium flash may never reach 
a state as the one given by our initial stellar model. Hence, this 
may influence the growth of the convection zone observed in our 
hydrodynamic simulations, as the thermodynamic conditions at 
the edges of the convection zone may differ. 

If a rapid growth of the convection zone indeed occurs, the 
main core helium flash studied here will never be followed by 
subsequent mini-flashes, as convection will lift the electron de- 
generacy in the helium core within 10 days. In addition, the he- 
lium core will likely experience an injection of hydrogen from 
the stellar envelope within a month and undergo a violent nu- 
clear burning phase powered by the CNO cycle. However, the 
growth of the convection zone within the core that is simulated 
in our models does not have to continue until it will reach the 
outer convection zone extending up to the surface of the star. 
Hence, mixing of nuclear ashes to the stellar atmosphere does 
not necessarily take place. But a fast dynamic growth of the in- 
ner convection zone will lead to a change of the composition of 
the stellar core (less carbon and oxygen), and consequently of 
the luminosities of low-mass stars on the horizontal branch. 

We found significant differences between the properties and 
the evolution of 2D and 3D models having a convection zone 
powered by (semi-)degenerate helium burning. However, as our 
3D models are likely not yet fully converged and as they cover 
only a relatively short period of evolutionary time (< 6000 s), 
long-term 3D simulations using a higher grid resolution are 
needed to obtain a better and more reliable understanding of the 
hydrodynamics of the core helium flash. 
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Appendix A: Hydrodynamic equations in spherical polar coordinates 

The hydrodynamic equations of a non-viscous multi-component reactive gas subject to a gravitational potential <P and having a heat 
conductivity K are given in spherical polar coordinates (r, 0, <p) by 

1 2 1 1 

d t p + -;d r (r pv r ) + — — dg(sm0pv e ) + — — <9<*,(pv ) = (A.l) 

1 o 2 1 1 Pvl P V l 

d,(pv r ) + -rd r (rpv r ) + — — -dg(smepv r Vg) + — — d (pv r v ) +d r p = -pd r <t> (A.2a) 

1 , , 1 , 1 pVaVr pvlcOS0 I p 

d,(pvg) + —d r (rpvf) + — — d g (sin6pv£) + ——d^(pv g v^) + — - + -d g p = — d 9 <$ (A.2b) 

1 o o 1 1 o PV<*V, pVdVgCOSO 1 p 

d t (p V<p ) + -d r (r 2 p V 2 r) + —^deisin OpvgVj,) + —Mpv\) + + H . + — -d*p = —3^ (A.2c) 

r l rsmQ rstn.8 * r rsinfl rsintf rsmd 

d,(pe) + -d r \r 2 [v r (pe + p) - K8 r T]} + -^—d B UnOMpe + p) - -d e T]\ + -^—d^ [v^pe + p) - -^—d^\ = (A.3) 

-p (v r 5 r <D + ^5 fl (D + -^3,*) + P£ (A.4) 

1,1 1 

d,(pX k ) + -rd r (r l pX k v r ) + — — <9 e (sin 6pX k v e ) + — — : d^(pX k v^ = pX k , k=l... N nuc (A.5) 

where p, v r , vg, v$, p, e, T, s, X k , and X k are the density, the radial velocity, the 0-velocity, the rotation velocity, the pressure, the 
total specific energy, the temperature, the energy generation rate per mass due to reactions, the mass fraction of species k, and the 
change of this mass fraction due to reactions, respectively. N nuc is the number of species the gas is composed of. 



Appendix B: Energy fluxes 



The various contributions to the total energy flux (Hurlburt et al. 1986; Achatz 1995) can be obtained by first integrating the 
hydrodynamic energy equation given in Appendix|A over the angular coordinates 6 and </>. Then, one decomposes both the specific 
enthalpy s + p/p (where s is the specific thermal energy) and the specific kinetic energy v,V;/2 into a horizontal mean and a 
perturbation, / = f + f , and obtains^] 

d l E + d r (F c + F K +F R +F E ) = (B.l) 



where 

E = | pedV 
Jv 



F K = 



v,p ■ |e+ — j r 2 dQ 



Kd,-T r 2 dQ 



i =1,2,3 



Fe — 4nr v, p ■ 



e + — + -VjVj + O 
p 2 



(B.2) 
(B.3) 
(B.4) 
(B.5) 
(B.6) 



Here, the gravitational potential <1> is assumed to be constant for simplicity. The sum of the various flux terms F, give the total 
energy transported per unit time across a sphere of radius r by different physical processes. One has the convective (or enthalpy) 
flux, Fc, the flux of kinetic energy, Fx, and the flux due to heat conduction and radiation, Fr. Finally, Fe, includes all terms causing 
a spherical mass flow, i.e., the model's expansion or contraction, while Fc and rest on deviations from this mean energy flow 
(vortices). The latter are the major contributors to the heat transport by convection. 

In a similar way one can also formulate a conservation equation for the mean horizontal kinetic energy that provides fur- 
ther insight into the effects of convective motions. Using the other hydrodynamic equations (Eqs. A.l to 
3f(pv,v,/2) = Vid t (pvi) - VjVidp/2, one finds 



A. 2c i, and the relation 



8 t E K + d r (F K + F P + F E w:) = P a + P p + P e .k 



(B.7) 



Note that the equations given in Mocak et al. ( 2008 ) that correspond to Eqs. i B.2 1, B.6 and i B.8 1 contain some small typographical errors 
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With Fk as introduced above, one obtains 



E K 
F P 
Fe,k 
Pa 

Pp 

Pe,k 



dV 



r 2 dQ 



!v2 ViV 

- (j) v r p' 

- ^) v r p 

47rr 2 ■ (pdm- v,.p5 r 0>) , i = 1,2,3 



~ " ~2 



ViV. 

+ 

P 

3,<D r 2 dQ 



(B.8) 

(B.9) 

(B.10) 

(B.ll) 

(B.12) 
(B.13) 



where the P, are source or sink terms of the kinetic energy. They are separated into the effect of buoyancy forces (Pa), and the 
work due to density fluctuations (Pp, volume changes). By analyzing the various P, one can determine what causes the braking or 
acceleration of the convective flow. The acoustic flux, F P , describes the vertical transport of density fluctuations. F E K and P E K 
describe the effect of expansion (volume work, and work against the gravitational potential), similar to F E in Eq. (B.6 1. 



